clear;

% ordering of firms:
% denmark: fringe, 5 Danes, Nordex
% germany: fringe, 5 Danes, 5 Germans 
%'distance_bonus','distance_nordtank','distance_micon','distance_vestas','distance_windworld',
    %'distance_nordex','distance_enercon','distance_fuhr','distance_suedwind','distance_tacke';


load data_structures
load est_point
%The project-to-project distance matrix
load dist_matrix_data

firmnames = {'Fringe', 'Bonus (DK)',  'Nordtank (DK)', 'Micon (DK)', 'Vestas (DK)', ...
    'WindWorld (DK)', 'Enercon (DE)', 'Fuhrlaender (DE)', 'Nordex (DE)', ...
    'Suedwind (DE)', 'Tacke (DE)'};

%Define Euler's Constant (not really needed, since its E(e), but just for
%completeness...
C = -psi(1);

%First seperate german and danish parts of the rho vector
rho_ger = e_rhohat(1: m.num_pro_ger*(m.num_firms_ger + 1),1);
rho_dnk = e_rhohat(m.num_pro_ger*(m.num_firms_ger + 1)+1 : m.num_pro_ger*(m.num_firms_ger + 1) + m.num_pro_dnk*(m.num_firms_dnk + 1),1);

%Now reshape rho so it is in firm x project form (like the Y matrix)
rho_ger = reshape(rho_ger,m.num_firms_ger + 1,m.num_pro_ger)'  ;  % num_pro_ger  x  (num_firms_ger + 1) matrix
rho_dnk = reshape(rho_dnk,m.num_firms_dnk + 1,m.num_pro_dnk)'  ;

%Get the generalized errors:
%  E[ek | y != k] = C - log(pk)  ; E[ek | y = k] = C + (pk/(1-pk))log(pk)
ge_ger = C - log(rho_ger).*(m.y_matrix_ger) + (rho_ger./(1-rho_ger)).*log(rho_ger).*(1-m.y_matrix_ger);
ge_dnk = C - log(rho_dnk).*(m.y_matrix_dnk) + (rho_dnk./(1-rho_dnk)).*log(rho_dnk).*(1-m.y_matrix_dnk);

%Here we need to define the spatial weight matrix across all projects:
nP = m.num_pro_ger + m.num_pro_dnk;
%nP = m.num_pro_dnk;

%%
% Define W:

for i=1:size(dist_matrix,1)
    for j=1:size(dist_matrix,1)
        if dist_matrix(i,j) == 0;
            dist_matrix(i,j) = 10;
        end
    end
end


%The "Correct" Distance Matrix, everything within a given distance
%Recall that the distance matrix is in km, 

%Euclidean Distance Cutoff
W = (dist_matrix < 30);
% Inverse Euclidean Distance...
%W = (1./ ((dist_matrix+1)) );
%W = (1./ ((dist_matrix)) );

% Inverse log Euclidean Distance...
%W = (1./log((1*dist_matrix)+2));

%This line just removes the diagonal.
W = W - diag(diag(W));

%Fake W:
%Used to confirm that the test is programmed correctly...if spatial weights
%are random there should be no spatial autocorellation.
%W = rand(nP);
%W = (W < .5);
%W = W - diag(diag(W));

%C pad generalized errors for firms that are not in denmark
ge_dnk = [ge_dnk C*ones(m.num_pro_dnk, m.num_firms_ger - m.num_firms_dnk)];
ge = [ge_ger ; ge_dnk];

%Let's start by checking the outside good:
tstRes = zeros(11,2);
out = zeros(11,3);
for f = 1:11

%Generealized residual demeaned, recall that E[e] = C.    
e = ge(:,f) - C;


%% This is the test...
moranI = e'*W*e;

Vm = 0;
for i = 1:nP
    for j = 1:nP
        Vm = Vm + ( (W(i,j)^2 + 2*W(i,j)*W(j,i))*e(i)^2*e(j)^2 );
    end
end
tstRes(f,1) = (moranI / sqrt(Vm));
tstRes(f,2) = 2*(1 - normcdf(tstRes(f,1)));
tstRes(f,3) = moranI / (e'*W'*W*e);
tstRes(f,4) = length(W)*(e'*W*e)*(sum(sum(W)))^(-1)*(sum((e-mean(e)).^2))^(-1);

%Now lets just use the heteroskedasticity robust regression test:
%This is the one we will report in the paper since it's the easiest
%to explain and directly gives us the small standard errors which
%we believe are driving the result.
X = [ones(size(e)) W*e];
beta = (X'*X)\(X'*e);

%This part calculates heteroskedasticity-robust White standard errors...
resid = e - X*beta;
B = zeros(2,2);
for nn = 1:length(e)
    B = B + X(nn,:)'*X(nn,:);
end
V = (X'*X)\ B /(X'*X);
se = diag(V).^.5;
out(f,:) = [beta(2) se(2) beta(2)/se(2)];

end

tstRes

%Use this to build the Appendix Table:
for f = 1:11
    fprintf('%s & %4.3f & %4.3f & %4.3f \\\\ \n', firmnames{f}, out(f,1), out(f,2), out(f,3));
end